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We study the numerical propagation of waves through future null infinity in a conformally com- 
pactified spacetime. We introduce an artificial cosmological constant, which allows us some control 
over the causal structure near null infinity. We exploit this freedom to ensure that all light cones 
are tilted outward in a region near null infinity, which allows us to impose excision-style boundary 
conditions in our finite difference code. In this preliminary study we consider electromagnetic waves 
propagating in a static, conformally compactified spacetime. 


I. INTRODUCTION 

The German noun “das All” refers to the Universe, 
but sometimes in a sense unfamiliar in English where 
it translates better as “outer space.” A spacecraft, af- 
ter launch, might be described as entering das All. In 
the computations reported here, it is precisely this outer 
space — the physical details far from a source of gravita- 
tional (or in our case electromagnetic) waves — which is 
excised from the computational model in analogy to the 
way the internal mysteries/singularities of a black hole 
are often excised in computational models. We develop 
and implement the proposals made earlier in Refs. [1-3] 
to obtain wave forms promptly at infinity by using hy- 
perboloidal time slices in a Cauchy evolution, to avoid 
infinitesimal time steps through conformal compactifica- 
tion, to convert outer boundary conditions to excision 
boundary conditions by using an artificial cosmological 
constant, and to defend the physical domain against com- 
putational boundary errors by tilting the light cones fur- 
ther outward in the region beyond the outer horizon. The 
aim of this paper is to explore the application of this set 
of boundary treatment tools in a low cost application in 
order to see whether there are unexpected impediments 
to its application, and to get guidance for the best ap- 
proaches for its application to larger problems. 

The use of hyperboloidal slicings in numerical relativ- 
ity has been carefully studied (e.g. Refs. [4-7]), and there 
are existing numerical results generated on compactified 
domains (e.g. Ref. [8]). This paper contributes to this 
literature by combining these two techniques in a space- 
time with an artificial cosmological structure designed to 
make it easier to construct boundary conditions at 
and by studying the effectiveness of modifying the light 
cone structure in the region of the analytically continued 
spacetime exterior to to prevent non-physical, in- 


coming radiation from piling up on the non-physical side 
of . That said, our approach to the problem does not 
require special modifications to the form of the differen- 
tial equations solved to model the physical problem. It 
should be compatible with exisiting numerical codes des- 
gined to solve the initial problem, described by some vari- 
ant of the the Einstein equations in an Arnowitt-Desner- 
Misner 3+1 decomposition [9]. 

Our work is also complementary to community wide 
efforts to use mesh refinement techniques in numerical 
relativity. Several groups have recently and successfully 
incorporated mesh refinement techniques into their Ein- 
stein solvers (i.e. Refs. [8, 10-18]). This technology allows 
groups to push their computational boundaries to very 
large coordinate distances, reducing the effects of incom- 
ing radiation and constraint violating modes introduced 
by applying boundary conditions at a finite distance. 

While mesh refinement has greatly increased the power 
of modern codes, it also has limitations. It is now fairly 
easy to surround a computational domain with coarser 
and coarser grids, but at some point the coarse grids 
cannot resolve a wavelength, see, e.g. Refs. [13] and [19, 
Section 4.3]. They can, however, insulate the active, wave 
generation regions of the computational domain from 
boundary problems for a limited period of time. During 
this period waveforms have been successfully extracted 
from inner regions as small as 2.4 wavelengths from the 
center [18]. The hyperboloidal slicings considered here 
and elsewhere have the advantage that they asymptote 
to the outgoing light cones in the spacetime. Outgoing 
waves, therefore, appear asymptotically constant and re- 
quire very few points to resolve them all the way to 

The rest of the paper is organized as follows: In Sec- 
tion II, we biefly summarize the formalism we use in our 
code, both for the background spacetime metric and for 
the Maxwell equations themselves. In Section III we de- 
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scribe the numerical methods used to solve the equations. 
We summarize the numerical results and discuss the im- 
plications for future work in Section IV. Additional detail 
is provided in the Appendix. 


II. FORMALISM 


A. Spacetime Metric 


As in Refs. [1, 2] the background metric for this study 
will be the de Sitter spacetime 

ds 2 =■ -dT 2 +dX 2 +dY 2 +dZ 2 + {R 2 /L 2 )(dT-dR) 2 (1) 

with an artificial cosmological constant A = 4-3 /L 2 . We 
are interested in the limit 1/L 2 — » 0 when this becomes 
Minkowski spacetime, and the cosmological constant is 
used to make boundary conditions numerically simpler 
(we hope) at the de Sitter horizon than they would be 
at flat spacetime’s The de Sitter horizon is the 

null hypersurface R — L in this metric. As before we 
introduce the coordinate changes 


T r 2 / 2 

— = u 4- r 777 

s 1 — r 2 / 4 

JC _ x l 

s 1 — r 2 / 4 


(2a) 

(2b) 


as used in Refs. [1-3], where R 2 = X 2 +Y 2 + Z 2 ~ X i X i 
and r 2 = x 2 +y 2 + z 2 = x l x\ This brings X + in to r = 2 
in the Minkowski case and in the de Sitter case makes 
this a spacelike hypersurface beyond the de Sitter hori- 
zon which is part of de Sitter spacetime’s future causal 
boundary which is also called The hypersurfaces of 
constant u are then hyperboloids 


(T~s(u-1)] 2 -R 2 = s 2 (3) 


in Minkowski spacetime and are also asymptotically null 
spacelike hypersurfaces in the de Sitter modification. 
This leads to a metric which is singular at r = 2, but only 
in a conformal factor s 2 /q 2 = s 2 /( 1 - r 2 / 4) 2 wilich does 
not appear in the Maxwell equations. Thus we can drop 
this conformal factor and our test problem is to solve the 
Maxwell equations in the resultant metric which is of the 
form 


ds 2 = - a 2 dt 2 4- 7 ij{dx l + 0 l dt)(dx j 4- 0 j dt ) (4) ' 


with t = u. When we choose W = 0, the following 
equations give the analytically continued and conformally 
regulated de Sitter metric described above: 


a 




Hj 


(i - WQ( i + j) 

\A +(f ) 2 


l- I y +(t ) 2 5 

X 9 

l+(f) 


fsS\ 2 

I J x i x j + Oij 


(5a) 


(5b) 


(5c) 


where 


_ (1 - W)r 1 

“ (1 + |) 2 + 4 



( 6 ) 


Figure 1 below shows the coordinate speeds for light rays 
on the inward side of the light cone v\ n , on the outward 
side of the light cone v out , and for a timelike center of 
the light cone u ce nter = -/? normal to the time slices of 
constant u. Because this allows light rays to move inward 



-beta (r) 
v_in (r) 


FIG. 1: Coordinate speeds for the de Sitter lightcone: inward 
(vin), normal to the spacelike time slices (—/?), and outward 
(v ou t). For the Minkowski metric the graph looks very similar, 
except that v\ n is there never positive. 


(^in < 0) in some of the region beyond but within 
our computational grid, we choose instead to take 


W(r) = 0(r - 2) 


(^~2) 4 

'(r-2) 2 + |] ; 


( 7 ) 


to be a smooth step function that modifies the metric 
beyond r — 2 in a way that keeps the “incoming” coor- 
dinate speed of light V{ n = dr/ du positive (that is outgo- 
ing) when r > 2. The function W is a C 3 modification 
of 0(r — 2), the Heaviside step function that jumps from 
0 to 1 at r = 2. Note that for r < 2, this metric has 
the surprising property that the outgoing speed of light 
v out = (1 -f-r/2) 2 is independent of the cosmological con- 
stant paramter L (cf. equation 8 in Ref. [1]), Although 
in general the presence of a cosmological constant alters 
the causal structure of a spacetime, and indeed it does 
modfiy the ingoing radial speed of light in the present 
example, it does not change the outgoing radial speed 
of light from that of the original Minkowski spacetime. 
This fact will prove useful for numerical diagnostics. Fig- 
ure 2 below shows the modified light cone directions in 
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— — — — -beta(r) 

— - v_in{r) 

FIG. 2: Coordinate speeds for the modified de Sitter (eqn. 7) 
iightcone: inward light rays, motion normal to the 

spacelike time slices (— j3) } and outward light rays (v ou t)- 


v_in (r) 



FIG. 3: Coordinate speeds inward (ui n ) for light rays in the 
modified (eqn. 7) de Sitter metric. 


this case. For this modified de Sitter metric using equa- 
tion (7) one finds very low inward light velocities in the 
region just outside as seen in Fig. 3. This can lead to 
fields propagating along the inward side of the Iightcone 
spending very long periods of time in the region r « 2.3 
as illustrated by the world line of a light ray in Fig. 4. 


dr/du = v_in {r ) 



FIG. 4: World line of a light ray on the inward side of the 
light cone in the modified de Sitter metric. 


B. Maxwell Equations 

We take the Maxwell equations in the form given in 
the DeserFest paper [1]. The quantities involved are all 
various components of the usual 4-dimensional Maxwell 
fields Fp V and = \^ 99 ^ a 9 u ^F Q p. In particular, for 
any 3 -F 1 metric of the form Eq. 4 we define, 

J 0 i = pi ; Fi = [ijk]H k (8a) 

and 

F i0 = Ei , F ij = [ijk]B k . (8b) 

Here [ijk] is the completely antisymmetric symbol [ijk] = 
0, ±1 with [123] = +1. The fundamental fields in the 
formulation are B x and V 1 . The constraints read 

diB 1 = 0 = BiV* . (9) 

These fields are each three dimensional vector densities 
so these equations are three dimensionally co variant. In 


terms of these primary fields, two auxiliary fields 
Hi are computed by the formulae 

Ei and 

Ei = + [ijk]0tB k 

(10a) 

and 


Hi = (a/V7)7 U& ~ [ijkWV k . 

(10b) 

The evolution equations are then 


do# = -[ijk -]djE k 

(Ha) 

and 


daTt = [ijk]diH k ■ 

(lib) 
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Note that the metric does not appear in equations (9) nor 
(11), and that the metric quantities which do appear in 
equations (10) are conformally invariant. This is a first 
order PDE system in flux conservation form. 

C. Initial Data 

For initial data and some testing purposes we have 
used the solution of the flat spacetime Maxwell equations 
employed by Knapp et al. [20] and subsequently by Fiske 
[21]. This solution is a wave pulse, which propagates 
smoothly from through the origin and out to . 
It consists of a torroidal V field and a poloidal B field 
generated from a vector potential, described in Ref. [1], 

A = AidX* = f sin 2 0d<f> (12) 

where 

sin 2 9d<t> = 4r ( XdY = YdX) = \ ( xdy - ydx) ' (13) 
R z r z 

and 

/ = (T - 2A U) e~ w2 - (T + 2Av) e" Ava . (14) 

In these equations U = T — R and V = T + R are the 
Minkowski retarded and advanced times, which are then 
expressed in terms of the compactified coordinates u and 
r as U/s = u-2 + 2/(l + ^r) and V/s = u-2-\-2/(l-~r). 
Finally, B x and V 1 are calculated by differentiating this 
A field. The term containing e~ xv can be smoothly 
omitted when r > 2 since all its derivatives vanish at 
r = 2, but the term with e~ xu must be retained to 
preserve the constraints, and we thus have nonzero initial 
data in the unphysical region r > 2, and did not have 
sufficient resolution to reduce the pulse width to make 
these unphysical data insignificant. 

For all of the simulations described in this paper, we 
choose A = 1 (or, equivalently, we measured distances in 
the uncompactified space in units of A -1 / 2 ). The pulse 
width is then A -1 / 2 . We also always choose s = 1 to set 
the size of the region near the origin where the hyper- 
boloidal slices are somewhat flat. Our results reported 
below used either s/L = 0 for the Minkowski case, or 
s/L = 0.1 for the de Sitter and modified (eqn. 7) de Sit- 
ter cases. Thus the r — const hypersurfaces were not 
spacelike for the de Sitter case when r > 2.21 

III. NUMERICS 

A key point in our approach to this problem is that we 
do not require special numerical techniques to handle the 
compactified spacetime at or within . At all interior 
points, we use second-order, finite-difference approxima- 
tions to all derivatives. At the computational boundary 
beyond we use a second-order, one-sided stencil for 



FIG. 5: The L2 norms (over the region r < 2) of the con- 
straints plotted vs. time. The higher resolution data is multi- 
plied by the factor appropriate to demonstrate second order 
convergence. The initial values arise from B and V fields 
which are exact analytic solutions of the constraints, eval- 
uated at the grid points and then differenced to form the 
constraints. 


normal derivatives. We integrate forward in time using 
the iterated Crank-Nicholson method commonly used in 
numerical relativity codes [22]. 

Our decision to use one-sided differencing at the com- 
putational boundary was (in our view) the simplest thing 
to do. It turns out, however, to be equivalent to the 
black hole excision boundary conditions first studied in 
Ref. [23]. To be precise, along a normal to the boundary, 
the composition of third order extrapolation to a cell just 
outside the boundary (as used for excision) with second- 
order center differencing at a cell just inside the boundary 
gives the same stencil as second-order one-sided differenc- 
ing at the cell just inside the boundary. This is consistent 
with our mental picture of excising a region exterior to 
and reflects the similarity between J 25 ' 4 ’ (at which all 
characteristics point outward) and an apparent horizon 
of a black hole (at which all characteristics point inward). 
As with black hole excision, it is critical for stable and 
accurate simulations that all light cones at the excision 
boundary tilt towards the excised region, hence the need 
for the beyond- J 25 ' 4 ' metric modification described in Sec- 
tion II A. Without these equation (7) modifications, the 
light cones would point outward only in a neighborhood 
of J^ 4- which did not extend out to the computational 
boundary of the cubic grid enclosing 

We have verified that our code converges at second or- 
der accuracy throughout the physical domain. We paid 
special attention to the region in the vicinity of J^ + , and 
we are satisfied that the code coverges there as well. Fig- 
ure 5 explicitly shows a two-point convergence test on 
the constraints over the physical region r < 2, which are 
seen to be second-order convergent to zero and nearly 
constant. 



0 
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FIG. 6: A space-time graph of r 2 B z . Notice that the wave 
passes through located at r = 2 in our coordinate system. 


IV. RESULTS 

When evolving with the original, conformal, de Sitter 
metric given in Eq. (1), large errors eventually appear 
near the computational grid boundary, grow exponen- 
tially in time, and propagate inside S rAr . This problem 
was slightly worse when the cosmological constant was 
omitted. These cataclysmic errors originated at the grid 
boundary which in these cases was not a spacelike hyper- 
surface for which excisions (at apparent horizons) were 
designed. [The (correctly zero) field B z took values near 
10 3 -i-ou/4 at e( jg es 0 f the z — 0 plane after u = 4, and 
values 10 s smaller at for 4 < u < 20.] As anticipated 
in Section II A we proceded to artificially tilt the light- 
cones outwards beyond to make the grid boundaries 
spacelike. 

Figure 6 shows a space-time view of the z component 
of the magnetic field for an evolution using the modified 
(eqn. 7) de Sitter metric. In order to make clear that the 
wave is cleanly passing through through the course of 
the evolution, we have multiplied the field by r 2 , making 
the scale of the field comparable at the orgin and at «/ + . 
This can be seen in a more local way in Figure 7, which 
shows the time evolution of B z at a single point on , 
for four cases. One is the analytic value of the field (in 
extended Minkowski spacetime) for comparison. Those 
which fail to keep the field zero after the pulse passes 
used the two (Minkowski and de Sitter) metrics where the 
equation (7) modifications to make the grid boundaries 
beyond spacelike were not implemented. 

Note that the analytic solution is with zero cosmolog- 
ical constant, whereas the numerically evolved solution 
is with finite cosmological constant, and yet the wave- 
forms coincide. This phenomenon can be understood as 
follows. Because, as pointed out in II A, the outgoing 
speed of light is the same in either case, the phase of 
outgoing waves can be said to be independent of a cos- 
mological constant here. Further, since by construction 
the initial wave pulse is identical in our Minkowski and 
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FIG. 7: The time evolution of B z at a point on j^ + . We show 
both a numerical solution generated by our code using the 
modified de Sitter (engineered lightcones beyond J?*) and the 
analytic solution with zero cosmological constant. The close 
agreement of the lines suggests that we are able to accurately 
capture the behavior of the fields as they pass through . 
Other numerical evolutions based on the original (compact- 
ified, analytically continued) Minkowski metric (s/L — 0.0) 
and the similiar de Sitter metric (s/L = 0.1). for which the 
grid boundaries were not spacelike, fail after t = u — 4. 


evolved with original metric, s/L=0.0 
evolved with original metric, s/L=0.1 
evolved with engineered metric. s/L=0.1 
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de Sitter simulations, by conservation of energy the wave 
amplitudes at are expected to be commensurate. 

These results show that many desired aims of this ap- 
proach are achieved. Waves propagate smoothly through 
and no difficulties appear to arise in the neighbor- 
hood of The compactified hyperboloidal slices allow 
waves to appear at a very short computational time 
after they originate in the central region. No modifi- 
cations were needed to integration algorithms designed 
for conventional spacelike Cauchy slicings. The excision 
at a cubic computational boundary, however, only main- 
tained acceptable behavior for a time several times the 
pulse width, and ultimately led to behaviors in the un- 
physical region beyond which became intolerabale. 
An example is given in Fig A2 where the constraint vio- 
lations beyond are seen eventually to increase expo- 
nentially. In addition, by using analytically constructed 
initial values satisfying the constraints, we had to accept 
nonzero initial conditions in the unphvsical region r > 2 
beyond , and the subsequent evolution of these fields 
complicated the computation beyond 

We believe that further developments of Cauchy evo- 
lutions with compactified hyperboloidal slicings and ar- 
tificial cosmology should be done with the computa- 
tional/excision boundary a spherical hypersurface at or 
modestly beyond The means for doing this have 

been developed, e.g. [24-27]. Then the possibility re- 
mains open that well posed wave equations could run 
indefinitely. As explained in the supplementary material 
(Appendix A), we did not here use the differential equa- 
tions from Section VI of [2] as our lightcone engineering 
(see Section II A) beyond makes the constraints there 
increase exponentially even from analytic arguments. A 
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spherical boundary could be spacelike at r > 2 without 
the need for such artificial tilting of the lightcones, and 
thus allow the use of the formulation [2. Equations 26] 
where constraints evolve causally. 
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APPENDIX A: AUXILIARY MATERIAL 

This section is mainly a collection of graphs which il- 
lustrate properties that are merely asserted in the main 
article. 


1. Constraints 

Even when the modified de Sitter metric is used, the 
solution in the unphysical region deteriorates as shown 
in Figs. Al and A2. 

In spite of these problems, the field at remains 
zero after the wave pulse passes (Fig. 7) for a moder- 
ately long time, which time can be lengthened by using 
higher resolution. See Fig. A3. Better control of the con- 
straints should be possible using the equation system 
[2, Equations 26] where the constraints propagate inside 
the lightcone and damp at a rate —difi 1 . However this is 
damping only for the Minkowski or de Sitter metrics as 
seen in Fig. A4 and becomes anti-damping for our modi- 
fication of the de Sitter metric (see Fig. A5 ) which makes 
our cubic grid boundaries spacelike. 

2. Initial conditions beyond 

A second reason which makes the use of the full cubic 
grid (rather than only a region inside a spherical bound- 
ary such as r = 2 or r - 2.1) inappropriate for numerical 
computation is the need for initial values beyond 


which is r = 2 or R = oo. 

With a spherical excision boundary and causal prop- 
agation of all modes, there is less room for bad initial 
conditions to intrude, and better reason to think that 
even poor initial conditions would be quickly flushed out 
of the computational domain. We have simply used as 
initial conditions the formulas for B and V as computed 
from equations (12) through (14) extended analytically 
by the coordinate transformation from TSYZ to uxyz 
of equations (2) . As a solution of the Maxwell equations 
in the Minkowski metric, these formulae give substan- 
tial activity in the r > 2 region, as seen in Fig.A6. In 
that Figure one solution has been modified by smoothly 
deleting the terms involving exp —XV 2 for r >2 which is 
possible since V/s = u — 2 + 2/(1 — ~r) becomes infinite 
at r = 2. However the pulse seen moving inward from 
beyond the grid will have depended upon initial condi- 
tions far outside r = 3.64 = 2.1\/3 which was the limit of 
our grid. Our initial data, however, are not evolved us- 
ing the Minkowski metric, but usually with the de Sitter 
or modified de Sitter metrics, where the decomposition 
into ingoing and outgoing parts of the initial conditions 
will differ seriously from the Minkowski case, especially 
well beyond . Thus we can understand that differ- 
ent, but equally visible, activity could occur in the r > 2 
region which would be actually solving the given differen- 
tial equations but be unrelated to the physical activity in 
the region r < 2. This appears to be occuring in the nu- 
merical results (modified de Sitter for which no analytic 
solution is available) shown in Fig. A7. 



FIG. Al: Samples of the time evolution of the constraint violation ]div£>! at z = 0 within our comptational cube \x l \ < 2.1 
whose edges are at r = 2.97 . 


r>2 



FIG. A2: The time evolution of an L2 norm of the constraint 
violations over the unphysical region r > 2 beyond within 
our computational cube |x z | < 2.1 whose corners are at r — 
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DivBeta (r) 



FIG. A5: The damping factor for constraints in the modi- 
fied de Sitter metric. Where div/? is positive, the constraints 
should grow exponentially at the plotted e-folding rate. 



FIG. A6: For the analytic solution in Minkowski spacetime: 
The magnetic field at 2 = 0 after the pulse from Fig 7 has 
exited the computational grid. 
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FIG. A7: A field component plotted along a line z = 0,x = 1 parallel to the y-axis. At the grid limit y = 2.1 one has r = 2.33, 
while y + is met at y = 1.73. Inside the field values appear to be convergent to zero after u = 4, while in the region beyond 
y = 1.85 or r = 2.1 a slowly moving pulse appears to remain independent of the selected resolution. 









